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ABSTRACT 

Aims. The first-order ordinary differential equation (ODE) that describes the mid-plane gravitational potential in flat 
finite size discs of surface density E(-R) oc R" (Hure & Hersant 2007) is solved exactly in terms of infinite series. 
Methods. The formal solution of the ODE is derived and then converted into a series representation by expanding the 
elliptic integral of the first kind over its modulus before analytical integration. 

Results. Inside the disc, the gravitational potential consists of three terms: a power law of radius R with index 1 + s, and 
two infinite series of the variables R and 1/R. The convergence of the series can be accelerated, enabling the construction 
of reliable approximations. At the lowest-order, the potential inside large astrophysical discs (s ~ — 1.5± 1) is described 
by a very simple formula whose accuracy (a few percent typically) is easily increased by considering successive orders 
through a recurrence. A basic algorithm is given. 

Conclusions. AppUcations concern all theoretical models and numerical simulations where the influence of disc gravity 
must be checked and/or reliably taken into account. 
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1. Introduction 

Gaseous discs in which the main physical quantities (den- 
sity, pressure, temperature, thickness, velocity) scale with 
cylindrical radius as power laws, i.e. "power-law discs" , 
represent an important class of theoretical systems. These 
are used customary to model accretion in ev olved bi- 
naries (IShakura fc Sunvae v 1973; Pringlc Il98lh . circum- 
stellar matter (jDubrullel 11992. : .Edgar. .2007h . the environ- 
ment of_jnassrve_UB;Ck_Jwles active galactic nu- 
cleiJConin-Souffrin fc DumontJ 119901 : iHurel 119981 : rSemera^ 
2004 ) or even the stellar component of some galaxies ( Evans! 



1994 IZhao et all Il999l) . In most applications however. 



power-law discs are truncated either to avoid diverging val- 
ues at the disc centre (such as density, mass) or in attempt- 
ing to reproduce the properties of observed discs of finite 
extension and mass. Although self-similarity is not compat- 
ible with the presence of edges, it is generally considered 
that power laws offer a good description of disc properties 
in some regions (far from the edges). Note that the presence 
of sharp edges c an be misleadin g when interpreting observa- 
tional data (e.g. lHughes et al.f2008i) . In general, the surface 
density S in the outer parts of discs is a decreasing function 
of the cylindrical radius R. Depending on the models, hy- 
potheses, aiid_objectSj_weJiavej_Jor instance, S oc i?^'^/^ in 
binaries (jShakura fc Sunvaev! [l973D . E oc i?±9/20 jn active 
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galactic nuclei (jCoUin-Souffrin fc Dumontlll990l ). S oc R~^ 
for a Mestel disc (lMestellll963D . or E oc i?-^/^ circum- 
stellar discs ()Pietu et al. l2007l ). In the context of station- 
ary viscous a-discs, a wide range of power-law exponents 
is allowed sinc e the temper ature T, E, and R satisfy the 
condition (e.g. lPringlelll98j] ): 

ETi?3/2 ^ ^g^^ 
while it is E oc i?"!/^ in /3-discs (|Hure et al.!l200l!l . 

The calculus of the gravitation potential of finite-size, 
power-law discs has received little attention yet. Several 
reasons can be put forward. Solving the Poisson equa- 
tion or computing the integral of the potential is not 
a trivial procedure, especially in the presence of edges. 
It is generally believed that gravity due to low mass 
discs is unimportant compared with that of a central 
proto-st ar or black hole, and c annot be probed (see 
however iBaruteau fc Massed 12008!) . Many studies employ 
the multi-pole expansion which is known to converge 
too slowly inside sources to be efficient for the numerica l 
applications (e.g . IClementI 11974 I Stone fc NormanI |1992[) . 
iHure fc Hersant! ()2007[ ) demonstrated that the mid-plane 
potential of flat power-law discs obeys an inhomogeneous 
first-order Ordinary Differential Equation (ODE). In this 
second paper, we discuss the exact solutions of this ODE 
for the entire physical range (outside and inside the disc) 
in terms of infinite series. In particular, it is shown that 
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the mid-plane potential is a combination of a power law 
for the radius R and two series of the variables R and 1 /R. 
Since these series converge rapidly inside large discs, it is 
possible to derive reliable approximations by truncating 
the series at low orders. 

This paper is organised as follows. The ODE for the po- 
tential is briefly recalled in Sect. [2] and its formal solution 
is derived in Sect. [3l In Sect. [H we express the potential 
at the two disc edges and consider a few special cases. The 
inside and outside solutions in the form of series are pre- 
sented in Sect. [S] In Sect. [6l we analyse the potential in 
the disc inside in detail, and in particular, the power-law 
contribution. Since all series involved converge rapidly, we 
are able to derive reliable approximations for the potential; 
this is done in Sect.[71 We discuss in Sect.[H]the case of discs 
with no inner and/or outer edge. The paper ends with a few 
concluding remarks. 
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Fig. 1. Configuration for a finite-size flat disc. 



2. The mid-plane potential in power-law discs from 
an Ordinary Differential Equation (ODE) 

Following [Hure k. HersantI (|2007[ ) (hereafter Paper I), the 
mid-plane potential tp due to a flat power-law disc satisfies 
the ODE: 



(1) 



where w = R/aout is the cylindrical radius in units of the 
radius Oout of the disc outer edge, s is the power-law in- 
dex of the surface density E, namely (generally, s < in 
astrophysical discs): 

S cx R\ (2) 

and S{'uj) is the piecewise defined function. Depending on 
the position in the disc (see Fig. [T|), we have: 

f.^[K(n7)-K(f)Ai+^] 

for < tu < A (region I), 



S{w) 



2Vout I n7K(n7) - K (^) 
for A < tn < 1 



(region II), 



(3) 



K(^)-K(i)A^+^ 

for w > 1 (region III), 

where A = a-m/aout < 1 is the axis ratio, Oin is the radius of 
the inner edge, -^out = 27rG'Eoutaout is a positive constant, 
Sout is the surface density at the disc outer edge, K is the 
complete elliptic integral of the first kind: 



K{x) = / 
Jo 



< a; < 1, 



(4) 



and G is the gravitation constant. The above ODE can 
in principle be solved in the entire radial domain since 
boundary conditions (both ip and the associated accelera- 
tion — droV') ^I'c known precisely at to = and at = oo for 
power law distributions. Although S is singular at the two 
disc edges (i.e. for w e {A, 1}), Eq. ^ is far more tractabl e 
in computing tp than the integral form fe.g. lDurandlll953l ): 



-4GI] 



out "out 



a + R 



-K 



2VaR 

a + R 



da, (5) 



whose integrand is logarithmically singular everywhere in- 
side the disc, or than the Poisson equation, which involves 
vertical gradients. 

3. Formal solution of the ODE 

A formal solution of Eq . ^ is found by setting (e.g. 
iRvbicki fc Lightman|[T979 



and 

S{zu) ~ vj 
The exact derivative of ^ is: 

d^ _ 
dw 



-(1 



-(1+^) 



V'(to), 

i+")5(n7). 



^-(1 + .)^' 
dw w 



(6) 
(7) 

(8) 



where we recognise, inside brackets, the function S. 
Therefore, we have: 



dtp ^ 

dvj 

whose formal solution for i]j is of the form: 

^{w) = '!/'(n7o) + / S{m')dw', 



(9) 



(10) 



Back-substituting ip and S from Eqs.® and we find 
the general expression for the mid-plane potential: 

Xn7o) , r S{m') 



dw' 



(11) 



^0 -J^a 

This solution is fully determined in the entire spatial do- 
main or part of it as soon as the potential is known 
at a given normalised radius wq. We observe that, for 
s 7^ —1, ^{w) is the mixture of a power law of the ra- 
dius (the first term in the right-hand-side) with exponent 
(e.g. lBisnovatvi-Koganlll975l:lEvans k. Readlll998D : 

p=l + s, (12) 

and a complicated function of the radius R (the definite 
integral) . 

In the following, we shall analyse Eq. pT|) analytically 
in terms of infinite series by conside ring in Eq. (|3t the ex- 
pansi on of K over its modulus (e.g. iGradshtevn fc RvzhikI 
I1965D : 

DO 

K(a;) = - X! ^'^^^ < 2; < 1, 

n=0 



70 = 1, 



(13) 



7«-i 



2n~l 
2n 



n > 1. 
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where 




1000 



n+1 



Fig. 2. Coefficient 7„ versus n. Terms a„ and \dn\ versus n 
for s = —1.5. 

4. Potential at the disc edges 

4.1. Inner edge 

To determine ij'i'cj) from Eq. (fTTjl for zu E [0, oo[, it is suffi- 
cient to calculate the potential at the two disc edgesQ, that 
is ip{A) and '(/'(I)- For this purpose, we use Eq. (O and 
define v — R/a < 1. After some algebr£0 we find: 



V;(A) = 4GSoutAi+^aout / - 

Ji V 



•t^.. (15) 



From Eq. (|T3|) . the potential at the inner edge is: 

oo 

^(A) = -Vout Ai+^ E ^" (1 - A2"-^-i) , (16) 

n=0 

where 



In 



2n - s - 1 



(17) 



and assuming A 7^ (see below). As 7„ ~ 1/n for large n 
(see Fig. [2]), consists of terms that vary asymptotically 
like ~ l/ri^. Since A < 1, ?/;(A) is a converging series. 

4.2. Outer edge 

At the outer edge, the potential is calculated in a similar 
manner, but using the variable u ~ a/R < 1. For H7 = 1, 
Eq. ^ writeJH 

V'(l) = -4GI]outaout / K{u)u^+'du. (18) 

J A 

Replacing K(m) by its series representation yields: 

00 

V-Cl) = -Vw E (1 - A2"+^+2) , (19) 



^ In contrast with the gravitational acceleration (the gradient 
of i/)), the potential is generally finite at the edges. 
^ In particular, we use the transformation: 



K^^^^ = {l+x)K{x), 0<x<l. 



(14) 



In 



(20) 



2n + s + 2 ' 

and by assuming A ^ (see below). As for and for 

the same reasons, "0(1) is also a converging series. 

4.3. Special values of s 

If the power law exponent s is such that: 

2n - s - 1 = (21) 

at a certain rank n = riA, then ?/'(A) must be, in practice, 
written in a slightly different form. This happens for s € £a 
with: 

= {-l,+l,+3,+5,...}. (22) 

Since 

1 - -r-9 

(23) 



1 - 2;9 

lim = — In a;. 



q^O q 

for a; > and any q, we have: 



lim /„(l-A2"-^-i) = -7„^lnA. 

n— >nA 



(24) 



Then, if s G S^, the potential at the disc inner edge is 
given by: 



V'(A) = -V-outA 



l+s 



J2 In (1 - A^-^-l) 



n=0 



(25) 



-jnA In A] . 

In a similar way, if s is such that: 

2n + s + 2 = (26) 

at a rank n = ni , then one term in Eq. (|19p must be treated 
separately. This happens for s E £1 where: 

fi ={-2,-4,-6,-8,...}. (27) 

From Eq. (|23p. we have: 



-0(1) = -Ipout 



E J„(l-A^"+^+^)-7nilnA 



(28) 



We note that s cannot belong simultaneously to set 5a 
and to set £1. 

4.4. Cases with A = 

The case A = occurs i) when the disc has no inner hole 
(i.e. fliii = 0) but finite size, and/or ii) when the disc has 
an inner edge but is infinitely extended (i.e. Cout — ^ 00). In 
the first case, iIj{A) (denoted -iAc in Paper I) becomes the 
potential at the origin of coordinates and has infinite value 
as soon as l + s < 0. In the second case, '0(1) represents the 
gravitational potential at infinity. It diverges if s -I- 1 > 0. 
Figure [3] summarises the ranges of s where the edge surface 
density, edge potential, and total disc mass are finite. We 
note that only discs having either i) Oin — 0, Oout ^ 00 with 
s > 0, or ii) flin > 0, Oout ^ 00 with s < —2 are physically 
meaningful since they are characterised by a finite surface 
density, a finite total mass and a finite pot ential. Mestel 
discs do not belon g to these categories (e.g. iMestell 119631: 
iHunter et aI]|1984D . 
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Fig. 3. Ranges of s where the edge surface density, edge 
potential and total disc mass are infinite. 

5. Solution of the ODE in the form of series 

We see from Eqs. ([3]) and (fT3l) that the function 5* can be 
easily expressed as a series. We have: 

in region I. 



■777^ ^ ^ 



(29) 



in region II, 
in region III. 

By inserting this general expression into Eq. we 
find for s ^ f A U f 1 and A ^ 0: 



+ V'out ^ 



n=0 



2n 



VJ 



l + s 



,l+s-2n 



^0 



(30) 



1 



-2n+l 



2n+2+s 
^0 



where the coefficients an and fe„ are respectively given by: 
1 — A^+''^^" in region I, 



1 in region II, 
in region III, 



(31) 



and 



hn — Jn ^ ^ 



in region I, 
^2ri+s+2 jj-j region II, 
^2n+s+2 _i region III. 



Since V'(^) ^nd ^/'(l) are available (see Sect. U]), we use 
vjQ = A and njp = 1 to simplify Eq. pO)) . Using Eqs. pB]) 
and (fT9l). we thus have: 



V'(tn) 

V'out 



l+s 



00 

E 

n=0 



^2«+l 



(33) 



total 

contribution of ttie power law 

0| , r 




Fig. 4. Potential {plain lines) in a power-law disc with axis 
ratio A = 0.1 and s = {—2.5, —1.5, —0.5}. The contribution 
of the power law, i.e. the term Aot^+^ in Eq. ([33|) . is shown 
in comparison. 



where: 



0, in region I, 







in region II, 



(34) 



with: 



0, in region III, 

-^n '-^n- 



(35) 

We note that A is a function of s. 

Figure [4] compares the total potential with the 

power law contribution (i.e. the term Aw^^'^) for three typ- 
ical values of the exponent s in a disc of axis ratio A — 0.1. 
We clearly see that, in a finite size disc: i) the gravitational 
potential is not a power-law function of the radius, ii) a 
power-law contribution is present inside the disc only, and 
iii) the power law is not the dominant part of the potential. 
As expected, spatial self-similarity is broken due to edges. 
We note that, outside the disc (i.e. in regions I and III), the 
series coincides with the multi-pole expansions. 

6. Potential inside the disc 

The determination of the gravitational potential is usu- 
ally straightforward outside the distribution wher e differ- 
ent k inds of expansions are efficient in practice (|Kellogd 
Il929[) . In contrast, it is problematic inside matter where 
the classical multi-pole approach fails to c onverge rapidly 
(e.g. IClementiHOTl IStone fc Normanlll992t) . 

6.1. Con verging series 

In region II, we have &„ = J„A^"+''+^. If we set: 



K VD 



(32) then 



-2n-|-l 



JnX 



2n-|-l^l-|-s 



< A 



l + s 



(36) 
(37) 



Figure [2] shows the term a„ versus n for s = —1.5 (in this 
case, /„ = J„). As a consequence, the three series involved 
in Eq. ([33]) converge rapidly since i) /„ and J„ both vary as 
1/n^ at large n, ii) all terms are positive for large n, and iii) 
zu < 1 and X < 1. This is interesting for the truncation of 
series and the construction of reliable approximations (see 
Sect. III). 
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Fig. 5. Coefficients A (open circles), A/^ (triangles), and 
Al (stars) versus the power-law exponent s of the sur- 
face density. A few remarkable exponents are shown 
(ope n squares): the Mestel (infinite) disc with s = — 1 
(|Mes tcl 1963), the H ayashi niqdel f or the solar neb- 
ula with s = —1.5 (|Havashil |1981[ ). and the a and 
/3-viscosity discs where s ft! —0.6 depe nding on mod- 
els ( Sliakura & Sun vaevI 119731: [Co Uin-Sou ffrin fc Dumont) 
Il99a lDubrullelll992tlRichard fc Zahnlll999[ ). 



6.2. The coefficient A 

The coefficient A is plotted in Fig. [5l It is symmetric with 

3 
2 



In 



4(4n-f 1) 



(4n-2s'-f-l)(4n + 2s' + l)' 



(38) 



s -I- |. For certain integer values of |s|, A is 



where s' 

strictly zero. This is in particular the case for s = — 3 and 
s = (see Appendix |A|) . As Eq. ([34| shows, A rises as soon 
as the exponent s is such that either 2n — s — 1 or 2n + s + 2 
is small (see Sect. |4]). Even, if n = (or n — ni), the 
coefficient A apparently contains a singular term, namely 
In^'cu^^^ (resp. Jni'cj^^^); however, this singularity exactly 
cancels with the term a„^n7^"^ (resp. fonjro^^"^"^). In 



practice, when s £ £a, Eq- 
we have: 



is no longer valid. Instead, 



■0out 



= (^A + 7riA ln^37) nj^^ 



a„n7 

ri=0 



C30 , 

n—0 



(39) 



where 



Aa = - 



oo 

E 

n=0 



oo 

/ -Vj 



(40) 



Similarly, if s e £i, the potential writes 
ip(w) 



(Al + 7„j In a:) w 

oc 



l-fs 



n=0 



E 



(41) 



2n+l ■ 



where 



E^« 



oo 

E 



Jn-i 



(42) 



A few values of A, , and Ai are listed in Appendix 
IbI We note that A a (or Ai) differs only from A by the term 
(resp. J„J. 

6.3. Convergence acceleration 

Once s is given, the coefficients A, A a, and Ai can be 
easily determined at the required accuracy. It is also pos- 
sible to improve the convergence rate of the associated se- 
ries. This accelerates the computation of the coefficients 
and makes their dependence with the exponent s more ex- 
plicit. Convergence acceleration is performed by using the 
properties of the definite integrals of the complete elliptic 
integrals of the first and second kinds. The demonstration 
reported in the Appendix [Cl vields. for s ^ £aU £i: 



A = 



1-6C 



- s(s + 2) ^ d„ - (s + l)(s + 3) ^ e„, (43) 



n=0 



n=0 



where 



4n2 - 1 ' 
4n2 - 1 ■ 



(44) 



and C is the Catalan constant (half the area under the 
function K). Numerically, the constant in A is 



1-6C 



-1.4310555380011220. 



(45) 



Figure [2] shows the coefficient c?„ for s = —1.5 (in this 
case, dn = e„). It follows that, for large n, d„ and e„ be- 
have like ~ 1/n'* asymptotically, and so A approaches its 
converged value far more rapidly than by means of Eq. (|34p . 

For s e £a, we then find: 



1-6C 
= 1- 



4nA 



--s(s + 2) dn 



n=0 



(46) 



(s + l)(,s + 3)^e. 



ri=0 



instead of Eq. (|40)l . and for s e £i, this is: 



1 — 6C 4ni -r-^ 

Al = + 7„, ,,^2 _ I - <^ + 2)}_^d„ 

n=0 



4n| 

(.s + l)(5 + 3) ^ e„. 



(47) 



instead of Eq. (|42l) . 

Depending on the exponent s of interest, a good approx- 
imation for A, Aa or Ai can be obtained by considering 
only the largest terms in the sum, i.e. all terms up to the 



rank rt « i(l-|-s) or n 



(s + 2). For astrophysical discs, 
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-1r 




-0.03 



-2 -1 
power law exponent s 



from Eq. (gg] 



s is around —1 meaning that we retain only the first term 
in Eq. (|43|) . We then find the following approximation: 



A: 



1-6C _ s{s + 2) (,9 + l)(s + 3) 

TT S+1 S + 2 



(48) 



whose accuracy is better than 3% for —2.7 < s < —0.3 as 
Fig.[6]shows. For s e {—2, —1}, we find from Eqs. (|46)) and 
(gZI): 

Aa = Ai^ (49) 

TT 

which is in good agreement with the converged value given 
in Table m 



7. Approximate formulae inside the disc 

Equations and (|¥T|) contain three rapidly converg- 

ing series that can be truncated to derive reliable approxi- 
mations for the potential. For s « — 1, only a few terms can 
be considered (see Sect. 15]). Although many truncations are 
possible, we have noticed that the most accurate approxi- 
mations of tp for disc^ are obtained provided the coefficient 
A (or Aa or Ai depending on s) takes its converged value. 
Under these circumstances, the A^-order approximation for 
the potential in region II becomes: 



N 



bn 



V'out 

which is, in its asymptotic limit 



N 



hm ^i^liw) ^ ^P{w). 



(50) 



(51) 



7.1. Zero-order approximation 



As argued in Sect. 16. 3i a reliable formula for the potential 
in astrophysical discs (for which s « —1.5 ± 1) is obtained 
by considering only the terms oq and bo/'cu, in addition to 



Fig. 6. Relative error in A when computed approximately 'S 




0,4 0,6 
normalized radius ns 

Fig. 7. The exact potential for a power law disc with 
s = —1.5 {thick line) compared to approximate values {thin 
lines) found from Eq. ([50]) for A^ = 0, 1, 2 and 3 (i.e. see 
Eqs. ([521), jM]), ([551), etc.). The largest errors are found 
around edges. 



the power law. At the lowest order, we thus havcQ: 
Ai+\nX 



V'ipp.(^) 

^out 



1 ifs = -2, 



= < 



Aa+Iuvu + X ifs = -l, 
1 



(52) 



Am'+' - 



s + 1 
72+^ 



otherwise. 



where, in this case, Ai — Aa ~ —1.386 (see Appendix [B|) . 
Figure [7] compares this zero-order approximation with the 
exact potential for typical disc parameters. It follows that 
the relative deviation |1 — tp^p./'>p\ does not exceed 10%, 
the deviation being the largest close to the edges. 

It is worth noting that the accuracy remains of the same 
order if coefficient A is determined by Eq. ([^5]) . This is con- 
venient when the explicit dependence of "0 on s is required. 



Here, discs are supposed to be objects of axis ratio A < 0.1. 



* A careful treatment of the singular cases s = { — 2, —1} shows 
that Eq. ([52};) yields Eqs. dUJi) and ([52b). 
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Under this hypothesis, the potential becomes (see note [3]): 

1 



ip{R) « ~27rGi;outaout x 



1 + s 



0.431 



1 



(l + 5)(2 + .) 



l + s 



1 



s + 2 \aout 




(53) 



Figure [5] shows the accuracy of this formula in the 
{w, s)— plane. We see that the relative deviation of ~ 10% 
observed previously for s = —1.5 holds globally for s 
roughly in the ranged [—3, 0]. This agrees with the fact that 
Eq. (I48p produces values of A within a few percents for this 
range of exponents. The deviation can be reduced at the in- 
ner and outer edges provided additional terms are included 
(see below). 

7.2. Higher orders 

If necessary, more accurate expressions are obtained by ac- 
counting gradually for following terms (each acting as a 
smaller and smaller correction). For iV = 1, we have: 



ipp.(^) 



rapp. 



■0out V'out 

and for N = 2, this is: 



+ 



4(1 -s) 4(4 + s) 



c 
o 
a. 

X 



o 
a. 




0.4 



0.6 



0.8 



normalized radius 



-2.5 



-1.5 



1.0 



-0.5 



(54) Fig. 8. Contour map showing the decimal logarithm of the 
relative error in the potential in the (ro, s)— plane when ip 
is approximated by Eq. ((53| . 



V'ipp.(^) 

V'out 



V'ipp.(^) 



64(3 -s) 64(6 + s)' 



(55) 



and so on. We note that Eq. ([55]) is particularly well suited 
to numerical computation since -0 can be determined by 
means of a recurrence procedure. A possible algorithm (not 
including the treatment of singular cases where s e S^USi) 
is proposed in the Appendix iDl 



8. Discs with no inner/outer edges 

8.1. Finite disc without inner hole 

If the disc has no inner edge but a finite size (i.e. oin = a nd 
Oout 7^ oo), then the ODE is (see iHure &: HersMrtl[200l : 



S{zu) 



'w'K{w), in region II, 



(56) 



K I 



in region III. 



It can be verified that the solution is still described by Eqs. 
((33l) and p3)). but the coefhcicnts a„ and 6„ are: 



and 



/„, in region II 
0, in region III 

0, in region II 
— J„, in region III 



(58) 



^ This range of exponents should be appropriate for most as- 
trophysical applications (see the introduction). 



8.2. Infinitely extended disc with inner hole 

If the disc is infinite but has an inner edge (i.e. ain > and 
Oout ^ oo), then 



where 



Y 



R vj 



is the new space variable, 



'7rF2 



'yK(y), in region I 
^K(^), in region II, 



(59) 
(60) 

(61) 



and i/jin = 27rGEinain is a constant (this is not the potential 
at the inner edge). The analogue of Eq. ([55]) is: 



''/'in ^ 



n=0 



y2ri+l 



(57) where A is still given by Eq. p4)) . 



and 



— /„, in region I 
0, in region II 

0, in region I 
J„, in region II. 



(62) 



(63) 



(64) 
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8.3. Infinitely extended disc 

If the disc is infinite, the ODE become homogeneous: 



0, 



and the solution is a power law: 



1+s 



(65) 



(66) 



where Rq is some reference radius. In this case only, a self- 
similar surface density can rigorously be associated with a 
self-similar potential. The presence of edges destroys this 
property. We note that, for s = —1 (i.e. Mestel's disc), the 
derivation of the ODE requires a careful treatment. The 
integral form, i.e. Eq. ([5]), gives: 



-4GSoao S 2C + - 



OO 

^ 2n 

,n=l 



1 



We then have 



d4 
dR 



1 

R' 



(67) 



(68) 



This expression is compatible with Eq. ([55k ) when flout — > 
OO (in this case, region III no longer exists and we have 

S = V'out/n7). 

9. Concluding remarks 

In this paper, we have derived an exact expression for the 
gravitational potential in the plane of flat power-law discs 
as a solution of the ODE reported in Paper I. This expres- 
sion is valid over the entire spatial domain and takes into 
account finite size effects. Inside the disc (the most difficult 
case to treat in general) , it consists of three terms of com- 
parable magnitude: a power law of the cylindrical radius R 
with index 1-1- s (where s is the exponent of the surface den- 
sity) and two series of R and 1/ R. In terms of convergence, 
our expression is by far superior to the multi-pole expan- 
sion method. Reliable approximations for the potential can 
be produced by performing fully controlled truncations. We 
have shown that the potential can be expressed by means 
of a simple function of R and s, which is valid to within a 
few percents in the range of exponents —3 ^ s < 0. This 
formula should be sufficiently accurate for most astrophys- 
ical applications. If necessary, more accurate formulae can 
be developped by including successive terms. These results 
should help in investigating various phenomena where disc 
gravity plays a significant role. 

An interesting point concerns the case of discs for which 
the surface density is not a power-law function. As shown 
in Paper I, it is easy to reproduce numerically the potential 
when the profile S(i?) is a mixture of power laws. From 
an analytical point of view however, the construction of a 
reliable formula for tp as compact as the one obtained here 
is not guaranteed at all. For instance, for an expansion of 
the form: 

M 



where s is an integer, each of the M + I series {fln,^n} 
should be truncated at a rank N > tia ^ {s + l)/2 (see 
Eqs. p3|) and ([50]) ). which corresponds to an approximate 
formula for i/j containing about 2M{M + 1) terms. The 
number of terms to consider can become prohibitively large 
when several power laws are required. 
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Appendix A: Vanishing coefficient A for s G {—3,0} 
For s E {—3,0}, the coefficient is : 



7n 



1 



2n + 2 2n - 1 



(A.l) 



We notice that the first sum is the definite integral of uK{u), which is 



E 



In 



2n + 2 



oo „i 

-^0 

„1 oo 

Jo 

/■I oo 

/ udu jnu'^^du, 

Jo 

— / 'K{u)udu, 

Jo 



(A.2) 



To fi nd the second sum, we com pute the definite integral of K{u)/v? in two ways. First, we have by direct integration 
(e.g. iGradshtevn fc Rvzhik|[T965[) : 



K{u) 



du = — 



TT ./n M 



E(u) 



E(l) - lim 



u^O U 



(A.3) 



2 


2 




h - hm 


TT 


TT ti^O 


2 






- lim — 


TT 


u— >0 U 



Second, by replacing above K by its series representation, we also find: 



2 KK 



du ■ 



TT ./n 



»1 oo 

oo 

^ 2n - 1 ^ 



= -70 



1 - hm - 

u^O u 



(A.4) 



70 



70 + 70 lim - + — + > - 

it-^o u \ — 1 2n — 1 

^ oo 

lim h — -^^^^ — . 

M^o u ^ 2n~l 



7o 



Since 70 = 1, we have 



E 7n _ _£ 
„ 2n- 1 ^ 



(A.5) 



and so A = for s = — 3 and s = 0. 
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Appendix B: Converged values of coefficient A, Ai, and for s G [—4,1] 



s 


-3 


— s 


s 






coefficient A 


relative error 


iteratioi 


-4 , 


, 00 


+1, 


,00 


-2. 


,50 


-9, 


. 65735902799760E-02 


+3 , 


. lE- 


■15 


6200 


-3 , 


, 95 


+0, 


.95 


-2. 


,45 


-5, 


. 08795971306598E+00 


+1 , 


. 6E- 


■16 


4763 


-3 , 


. 90 


+0, 


.90 


-2. 


,40 


-2, 


. 57897148688964E+00 


+9 , 


. 9E- 


■17 


6289 


-3, 


.85 


+0, 


.85 


-2. 


,35 


-1, 


. 73622548292480E+00 


+3 


. 9E- 


■17 


8653 


-3 , 


. 80 


+0, 


.80 


-2. 


,30 


-1, 


. 30966680644896E+00 


+2 


. lE- 


■17 


10760 


-3 , 


. 75 


+0, 


.75 


-2. 


,25 


-1, 


. 04923490407066E+00 


+8 


. lE- 


■17 


7977 


-3 , 


. 70 


+0, 


.70 


-2. 


,20 


-8, 


.71529232712215E-01 


+1 


. 4E- 


■16 


7154 


-3 , 


. 65 


+0, 


.65 


-2. 


, 15 


-7, 


.40760471984439E-01 


+2 


. lE- 


■16 


6696 


-3 , 


. 60 


+0, 


.60 


-2. 


,10 


-6, 


.38987239513759E-01 


+2 


. 8E- 


■16 


6379 


-3 , 


. 55 


+0, 


.55 


-2. 


,05 


-5, 


. 56193814181879E-01 


+3 


. 5E- 


■16 


6133 


-3 , 


. 50 


+0, 


.50 


-2. 


,00 


-4, 


.86319914494768E-01 


+4, 


. 3E- 


■16 


5926 


-3 , 


.45 


+0, 


.45 


-1, 


,95 


-4, 


. 25454571886386E-01 


+5 


. 2E- 


■16 


5743 


-3 , 


.40 


+0, 


.40 


-1, 


,90 


-3, 


. 70931368062394E-01 


+6 


. 3E- 


■16 


5576 


-3 , 


. 35 


+0, 


.35 


-1. 


,85 


-3, 


. 20839275213352E-01 


+7 


. 7E- 


■16 


5420 


-3 , 


. 30 


+0, 


.30 


-1, 


,80 


-2, 


. 73740796314985E-01 


+9 


. 4E- 


■16 


5270 


-3 , 


. 25 


+0, 


.25 


-1, 


,75 


-2, 


. 28500045775376E-01 


+1 


. 2E- 


■15 


5124 


-3 , 


. 20 


+0, 


.20 


-1. 


,70 


-1, 


. 84171894656878E-01 


+1 


. 5E- 


■15 


4980 


-3 , 


. 15 


+0, 


. 15 


-1. 


,65 


-1, 


. 39926048201652E-01 


+2 , 


. OE- 


■15 


4837 


-3 , 


. 10 


+0, 


.10 


-1, 


,60 


-9, 


.49911954664950E-02 


+3, 


. lE- 


■15 


4694 


-3 , 


. 05 


+0, 


.05 


-1, 


,55 


-4, 


. 86101640289562E-02 


+6 , 


. 3E- 


■15 


4550 


-3 , 


. GO 


+0, 


.00 


-1. 


,50 


-8, 


. 70956628336228E-15 


+3 


. 6E- 


■02 


4403 


-2 , 


. 95 


-0, 


.05 


-1. 


,45 


+5, 


. 16876888183838E-02 


+6 


. 4E- 


■15 


4252 


-2 , 


. 90 


-0, 


.10 


-1, 


,40 


+1, 


.07410381098617E-01 


+3 


. 2E- 


■15 


4098 


-2 , 


. 85 


-0, 


. 15 


-1, 


,35 


+1, 


. 68286941322651E-01 


+2 


. lE- 


■15 


3937 


-2 , 


. 80 


-0, 


.20 


-1. 


,30 


+2, 


.35664501680007E-01 


+1 


. 6E- 


■15 


3771 


-2 , 


. 75 


-0, 


.25 


-1. 


,25 


+3, 


. 11208301380443E-01 


+1 


. 2E- 


■15 


3596 


-2 , 


. 70 


-0, 


.30 


-1, 


,20 


+3, 


.97027142932255E-01 


+1 


. OE- 


■15 


3410 


-2 , 


. 65 


-0, 


.35 


-1, 


, 15 


+4, 


. 95854560428368E-01 


+8, 


. 6E- 


■16 


3212 


-2 , 


. 60 


-0, 


.40 


-1. 


,10 


+6, 


. 11318936974807E-01 


+7, 


.4E- 


■16 


2997 


-2 , 


. 55 


-0, 


.45 


-1. 


,05 


+7, 


. 48359798181055E-01 


+6 


. 5E- 


■16 


2759 


-2 , 


. 50 


-0, 


.50 


-1. 


,00 


+9, 


. 13893162088934E-01 


+5 , 


. 7E- 


■16 


2486 


-2 , 


.45 


-0, 


.55 


-0. 


,95 


+1, 


. 11792019234452E+00 


+5 


. lE- 


■16 


2151 


-2 , 


.40 


-0, 


.60 


-0. 


,90 


+1, 


. 37546760819327E+00 


+4 


. 5E- 


■16 


1662 


-2 , 


. 35 


-0, 


.65 


-0. 


,85 


+1, 


. 71019225943660E+00 


+4 


. lE- 


■16 


1335 


-2 , 


. 30 


-0, 


.70 


-0. 


,80 


+2, 


. 16159222263118E+00 


+3 


. 7E- 


■16 


1821 


-2 , 


. 25 


-0, 


.75 


-0. 


,75 


+2, 


. 80087471242395E+00 


+3 


. 4E- 


■16 


1997 


-2 , 


. 20 


-0, 


.80 


-0. 


,70 


+3, 


.77063202688011E+00 


+3 


. lE- 


■16 


2058 


-2 , 


. 15 


-0, 


.85 


-0. 


,65 


+5, 


. 40388067136072E+00 


+2 


. 8E- 


■16 


2036 


-2 , 


. 10 


-0, 


.90 


-0. 


,60 


+8, 


. 70024086898869E+00 


+2 


. 6E- 


■16 


1928 


-2 , 


. 05 


-0, 


.95 


-0. 


,55 


+1, 


. 86592556634413E+01 


+2 


. 4E- 


■16 


1681 


-2, 


.00 


-1, 


.00 


-0. 


,50 


-1, 


.38629436111989E+00 


+7, 


. 2E- 


■18 


7955 




95 


-1, 


.05 


-0, 


,45 


-2 


. 14370823367682E+01 


+2 


lE- 


16 


1770 


-1, 


.90 


-1, 


. 10 


-0, 


,40 


-1 


. 14939346514465E+01 


+1 


.9E- 


■16 


2141 


-1, 


.85 


-1, 


.15 


-0, 


,35 


-8 


.22454216651923E+00 


+1 


.8E- 


■16 


2399 


-1, 


.80 


-1, 


.20 


-0, 


,30 


-6 


.63018820764787E+00 


+1 


.7E- 


■16 


2598 


-1, 


.75 


-1, 


.25 


-0, 


,25 


-5 


.71250114438470E+00 


+1 


.7E- 


■16 


2755 


-1, 


.70 


-1, 


.30 


-0, 


,20 


-5 


. 14024384191494E+00 


+1 


.6E- 


■16 


2878 


-1, 


.65 


-1, 


.35 


-0, 


,15 


-4 


.77330268633758E+00 


+1 


.6E- 


■16 


2972 


-1, 


.60 


-1, 


.40 


-0, 


,10 


-4 


.54390925876444E+00 


+1 


.5E- 


■16 


3039 


-1, 


.55 


-1, 


.45 


-0, 


,05 


-4 


.41737401180816E+00 


+1 


.5E- 


■16 


3078 


-1, 


.50 


-1, 


.50 


+0, 


,00 


-4 


.37687923045294E+00 


+1 


.5E- 


■16 


3091 


-1, 


.45 


-1, 


.55 


+0. 


,05 


-4 


.41737401180816E+00 


+1 


.5E- 


■16 


3078 



comment 



values of Aa and Ai 



strictly (see Appendix 



values of Aa and Ai 



Table B.l. Values of A (4th column) computed within the computer precision (double precision; about 16-digit) for 
different power law exponents s or s' (columns 1 to 3) including those relevant for astrophysical applications (see also 
Fig. [5]). Also given are the relative error on the coefficient (5-th column) and the number of iterations (6-th column) 
required. 
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Appendix C: Convergence acceleration for A 

The coefBcient A i s Riven by a series whose con vergence can be accelerated by considering an interesting property of the 
integral of K fe.g. iGradshtevn &: Rvzhi3ll965D . namely : 



10 

where C is Catalan's constant. We can then write 

In 



„1 oo 

Jo ^ ' 2 4- 2n + 1 



(C.l) 



^ 1 + s - 2n \1 + s - 2n 2n + 1 J ^ 2n + 1 

n—0 n—0 



Gs + 2)^ 



-2n 2n 

In 



n=0 

AC 



^ (2n + l)(l + s-2n) tt 

Tl— 



(C.2) 



where terms vary like for large n. A second convergence acceleration can be obtained by considering the complete 
elliptic integral of the second kind : 



E(a;) = 

whose definite integral over the modulus x is 



7r/2 



l-x'^sm^A, 0<x<l, 



/ E(x)dx = --y- — ^ — - = - 

Jo 2^ (2n+l)(2n-l) 2 



C. 



We then have 



E 



In 



1 



1 



(2n+l)(l + s-2n) (2?! + l)(2n - 1) 
7" ,2/1 



E 



^^^(2n + l)(2n-l) 



(2n + l)(l + s-2n)(2n-l) n \2 



ra=0 



C 



It follows that : 



^ 1 + s - 2n ^ ' ^ 



In 



2 + s + 2sC 



n=0 



(2n + l)(l + s-2n)(2n-l) 



The second term in Eq. gives : 

oo 



2n + s + 2 

n—0 n—0 



1 



-G^ + i)E 



2n + s + 2 2n 

In 



1 \ °° 

+ 1 ^ 2n+l 



n=0 



AC 



(2n + l)(2n + s + 2) tt 



But, from Eq. (jC.4p . we have : 



E 



(2n+l)(2n + s + 2) 

n—0 ^ ^ ^ ^ n—0 



E' 



1 



1 



(2n+ l)(2n + s + 2) (2n + l)(2n - 1) 
In 2/1 



(2n+ l)(2n- l)(2n + s + 2) n \2 



7n 

^^(2n+l)(2n-l) 



C 



and then : 

V , =(g + l)(. + 3)V 

Finally, we have : 



In 



l + .s + 2(s + 3)C 



2n + s + 2 '"^ ' ' ^ (2n+ l)(2n- l)(2n + s + 2) 

n—0 ^ ^ ^ ^ ^ ^ 



(C.3) 



(C.4) 



(C.5) 



(C.6) 



(C.7) 



(C.8) 



(C.9) 



A^^Cn 

oo 

= s{s + 2)y 



n=0 



{2n + l)(l + s-2n)(2n- 1) 



— -Gs + l)(,s + 3)^ 



In 



1 -6C 



(2n+ l)(2n- l)(2n + s + 2) tt 



(C.IO) 
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Appendix D: A basic algorithm 

A possible algorithm is the following (see the text for cases where s E £i U £a)- 
1. disc parameters (see Sect. [5]): 

- set the inner edge Oin =??? 

- set the outer edge Uout =??? 



set the outer surface density E, 



out 



-777 



- set the power law exponent s =??? 
— *■ compute A 

2. power-law coefficient A : 

— *■ compute A (or Aa or Ai depending on s) at the computer accuracy (see Sect. [7]); see Eqs. (|43| . (|46|) and 
(|47p. or use pre- computed converged values (see Appendix [B]) 

3. radius 

- set the radius R =111 

— > compute dimensionless variables vo and X (depending on the region I, II or III) 

4. initializations 

- set 7o to 1; see Eqs. ([T^ 

- set ao to see Eq. (I5T]) 

- set 6o to 2^ ; see Eq. ([M]) 

- set -0 to 

5. main loop on n 

— > compute 7„ from 7„_i 
— > compute a„ from a„_i 
— > compute hn from 
^ compute nj^" and 

^ update V to V' + fln^^^" + &,in7^^"~^; see Eq. ([50)1 

6. power-law contribution 

^ update V to V' + v4ct7^+^; see Eq. ([50)) 

7. final potential value 

change ior ij] x 27rG'Soutaout 
The loop ends after N steps; the accuracy of the potential value is then given by the next term (rank N + 1). 



